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Abstract 

We present two methods for solving the electrostatics of point charges and multipoles on the sur¬ 
face of a sphere, i.e. in the space S 2 , with applications to numerical simulations of two-dimensional 
polar fluids. 

In the first approach, point charges are associated with uniform neutralizing backgrounds to 
form neutral pseudo-charges, while, in the second, one instead considers bi-charges, i.e. dumbells 
of antipodal point charges of opposite signs. We establish the expressions of the electric potentials of 
pseudo- and bi-charges as isotropic solutions of the Laplace-Beltrami equation in 1 S 2 . A multipolar 
expansion of pseudo- and bi-charge potentials leads to the electric potentials of mono- and bi¬ 
multipoles respectively. These potentials constitute non-isotropic solutions of the Laplace-Beltrami 
equation the general solution of which in spherical coordinates is recast under a new appealing form. 

We then focus on the case of mono- and bi-dipoles and build the theory of dielectric media in 
S 2 ■ We notably obtain the expression of the static dielectric constant of a uniform isotropic polar 
fluid living in £2 in term of the polarization fluctuations of subdomains of £ 2 - We also derive 
the long range behavior of the equilibrium pair correlation function under the assumption that it 
is governed by macroscopic electrostatics. These theoretical developments find their application 
in Monte Carlo simulations of the 2D fluid of dipolar hard spheres. Some preliminary numerical 
experiments are discussed with a special emphasis on finite size effects, a careful study of the 
thermodynamic limit, and a check of the theoretical predictions for the asymptotic behavior of the 
pair correlation function. 
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I. INTRODUCTION 


The idea of using the two dimensional (2D) surface of a sphere, i.e. the space S 2 , to 
perform numerical simulations of a 2D fluid phase can be tracked back to a paper by J.-P. 
Hansen et al. devoted to a study of the electron gas at the surface of liquid Helium [1], 
Subsequently, the same idea was used to study the crystallization of the 2D one-component 
plasma (with oc logr interactions) j^], to establish the phase diagram of the 2D Coulomb 
gas jij, and to determine some thermodynamic and structural properties of the 2D polar 
fluid (with oc 1/r 2 interactions) in its liquid phase ( 4 , 5 ]. The generalizations to 3D systems, 
implying the use of the surface of a 4 D hypersphere, i.e. the space S 3 , is due to Caillol and 
Levesque js]. Some improvements on this early work, mostly applications concerning 3D 
polar fluids, were published recently Jt|. 

Several simple ideas can be put forward to justify the use of a hyperspherical geometry 
in numerical simulations of plasmas and Coulombic fluids (i.e. fluids made of ions and/or 
dipoles) : 

• (i) The n-dimensional non-Euclidian space S n (in practice n = 2,3), albeit finite, is 
homogeneous and isotropic, in the sense that it is invariant under the group 0 (n-\-l) 
of the (n + 1 )D rotations of the Euclidian space E n+l ; it is thus well suited for the 
simulation of a fluid phase (in the bulk). 

• (ii) The laws of electrostatics can easily be obtained in S n and the Green’s function of 
Laplace-Beltrami equation (i.e. Coulomb potential) is known It has a very simple 
analytical expression, tailor-made for numerical evaluations. 

• (iii) The inclusion of charged or uncharged walls can easily be done in order to study 


the structure of liquids, plasmas, or colloids at interfaces 


12j 




The second point (ii) was partly overlooked by the authors of Ref. |4j, l5| who used an em¬ 
pirical 2D dipole/dipole interaction, unfortunately not deduced from a solution of Laplace- 
Beltrami equation in So, which casts some doubts on the validity of their results concerning 
the 2D Stockmayer fluid, notably those concerning its dielectric properties. 

Some years after we are now in position to propose in this paper two possible, both 
rigorously correct, dipole/dipole pair-potentials and to use both of them into actual numer¬ 
ical simulations of a 2 D dipolar hard sphere (DHS) fluid. Here, not only we extend to <S 2 
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the recent developments on 3D polar fluids in S 3 discussed in Ref. |]?j], but we also present 
additional new results on multipolar expansions in S 2 , not yet available in S 3 , in Sec. Q4] 
This article is thus a contribution to the physics of 2D fluids, colloids or plasmas, the 
atoms or molecules of which interact via electrostatic pair potentials derived from a solution 
of the 2D Laplace equation. We emphasize that the systems that we consider cannot be seen 
as thin layers of real 3D systems. In this case, the electrostatic interactions should be derived 
from the solutions of the 3D Laplace equation. Simulations of such 3D systems constrained 
to live in a 2D geometry are nevertheless also possible in spherical geometries. In that case 
one should consider a 3D fluid of S 3 made of particles interacting by pair potentials deduced 
from solutions of Laplace-Beltrami equation in S 3 and constrained to stay on the equator of 
the hypersphere. This geometrical locus indeed identifies with the sphere S 2 - This approach 
was used for the 3D restricted primitive model of electrolytes in Ref. 13]. 


In this work our aim is to demonstate the validity of the sphere S 2 to perform Monte 
Carlo (MC) simulations of a 2D polar fluid. We note that, although 2D dipolar fluids do 
not exist per se in nature, the model could be used via various mappings for applications 


as, recently, for the hydrodynamics of two-dimensional microfluids of droplets 14]. 


Our paper is organized as follows. After this introduction (Sec. ID) we discuss in depth 
the electrostatics of the distributions of charges in S 2 and their multipole expansions. The 
material reported in Sec. |TT] is an intricate mixture of old and new stuff. The underlying idea 
is the remark of Landau and Lifchitz in Ref. 15] that, in a finite space such as S 2 , Laplace- 
Beltrami equation admits solutions if and only if the total electric charge of the space is 
equal to zero. Therefore, in S 2 , the building brick of electrostatics cannot be a single point 
charge o as in the ordinary 2D Euclidian space E 2 . We are led instead to consider a pseudo¬ 
charge |6j, a neologism denoting the association of a point charge and a uniform neutralizing 
background of opposite charge. Alternatively we can consider a bi-charge, i.e. a neutral 
dumbell made of two antipodal charges of opposite signs +q and — q as in Ref. 16]. Both 
approaches yield independent electrostatics theories which are sketched in Sec. ITTl and UTT1 
The case of dipoles is then examined more specifically in Sec. IIII1 We need make a 
distinction between mono-dipoles, which are obtained as the leading term of a multipole 
expansion of a neutral set of pseudo-charges, and bi-dipoles, obtained in a similar way from 
bi-charges. As a consequence two possible microscopic models of dielectric media can be 
worked out, those made of mono-dipoles living on the entire sphere S 2 and those made of 
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bi-dipoles living on the northern hemisphere S^- 

These two types of models are then studied in the framework of Fulton’s theory 


17 
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in Sec. m Fulton’s theory realizes a harmonious synthesis of linear response theory and the 
macroscopic theory of dielectric media. It allows us to obtain 

• a family of formula for the dielectric constant e of the fluid, well adapted for its 
determination in MC simulations. 

• an expression for the asymptotic pair correlation in <S 2 under the assumption that it 
should be dictated by the laws of macroscopic electrostatics. 

Monte-Carlo (MC) simulations of the 2D DHS fluid in the isotropic fluid phase are then 
discussed in Sec. El We have chosen to report MC data only for two thermodynamic states, 
both in the isotropic fluid phase and extensive MC simulations of the 2D DHS fluid will 
be published elsewhere. These states could serve as benchmarks for future numeric studies. 


21| by means of MC simulations 


One of these states was studied many years ago in Ref. 
within standard periodic boundary conditions (see also (22j). Comparisons of our results 
with those obtained in this pioneer work are OK. A carefull finite size scaling study of our 
data yields results of a high precision, notably for the energy, probably unattainable by 
standard simulation methods. 

Conclusions are drawn in Sec (IVIf) . 


II. ELECTROSTATICS OF 2D CHARGES AND MULTIPOLES 


A. The Plane 


Let us first recall that the electrostatic potential Ve 2 (p) at a source-free point p = (x,y) 
of the Euclidian plane E 2 satisfies the 2D Laplace equation Ae 2 Ve 2 = 0, which in polar 
coordinates (p, p) reads 


Ae 2 Ve 2 — 


ld_ 

pdp 


dV\ 1 d 2 V 

^ dp J p 2 dp 2 


The general solution of Eq. (JT|) is 23, 24] 


( 1 ) 


oo oo 

Ve 2 (p, p) = a 0 + b 0 log p + ^ a n p n cos (np + a n ) + ^ b n p~ n cos (np + /3 n ) , (2) 

n =1 n =1 


4 












where a n , b n , a n , and f3 n are arbitrary constants. The terms of the r.h.s. of (J2]) which are 
singular at p = 0 may be interpreted as the potentials created by point multipoles located 
at the origin O, while those singular at p = oo as the potentials of point multipoles at 
infinity. In E 2 the potential of a point charge q located at O is — q log p up to an additional 
constant; it is the Green’s function of Eq. ([!]) in E 2 (without boundaries) and it satisfies the 
2D Poisson’s equation 


Ae 2 ( log p) = -2 tt6e 2 (M,M 0 ) = -2tcS( x)5(y) , (3) 


where S(x) denotes the ID Dirac’s distribution. 

Finally, the Green’s function — log(|p — p |) can be expanded in polar coordinates as [231 

"logflp-p'D = -logp< cos (n {^P — V 7 ) J j (4) 

where p< = inf(p, p) and p> = sup(p, p). 

Eq. (HD serves as a starting point for the 2D multipolar expansion of Ref. 24], Let us 
consider N charges qt of polar coordinates (pi, pi) and a point M (coordinates (p, p)) of the 
plane E 2 such that p > pi, V*. Making use of (HD one finds 

N 

Ve 2 (M) = -J^g*log(|p»pi|) 


2=1 


2 _, °°' _^ ^ 

V ° ~~ 9 logp + 9 —Q™ exp (-ini/ip) , 


2 ■ t —' ' np" 

71=1 V 


( 5 ) 


where Vo is some unessential constant. In (J5]) the index v assumes the values —1 and +1 
only. The complex multipolar moments are defined as 

N 

Qnv = exp (invipi) . (6) 

1=1 

Since Q n -u = Q* n u there are at most two independent multipole components at a given 
order n. The potential of a 2D multipole of order n therefore decays as 1 /p n . This remark 
justifies a posteriori the interpretation that we gave for the terms of the r.h.s. of Eq. (J2J). 


B. The Sphere 

Here, we examine how the basic laws of 2D Euclidian electrostatics swiftly evoked in 
Sec. Ill Al are changed when the system is wrapped on the surface of a sphere. Let us 
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first introduce some notations and recall some elementary mathematics. We denote by 
d>2(0, R) the sphere of center O and radius R of the usual 3D geometry, ft is a (Riemannian) 
manifold of the 3 D Euclidean space E 3 that we identify with M 3 . When we deal with the 
sphere of unit radius we adopt the uncluttered notation S 2 = S 2 (0,R = 1). Let M be the 
running point of S 2 (0,R), we define the unit vector z G S 2 as OM = Rz. The spherical 
coordinates are defined as usual : z = (sin# costp, sin6* sin <p, cos9) T with 0 < 0 < tt and 
0 < ip < 2n. The differential vector dz = d 6 eg + sin Odcpe^ allows us to obtain two 
orthogonal unit vectors (eg, e^,) that span the plane T 2 (M) tangent to S 2 (0, R ) at point M. 
Recall that e e = (cos# cos<p, cos 6 sin <p, — sin#) T and = (— sin <^ 3 , cos 99 , 0) T . In addition 
and quite specifically since this notion cannot be generalized to higher dimensions, one has 
e v = z x eg where the symbol x denotes the 3 D vectorial product. Finally the infinitesimal 
surface element is dS = R 2 dVt where the infinitesimal solid angle dVt = sin QdQdip in spherical 
coordinates. 

In the sequel we will make a repeated use of the unit dyadic tensor Us 2 (z) = egeg + 
of the plane T 2 (M). Note that the unit dyadic tensor of Euclidian space E 3 is given by 
Ue 3 = U,s 2 (z) + zz. It is a constant tensor independent of point M. These admittedly 
old-fashioned objects however allow an easy definition of the gradient in S 2 (0,R), or first 
differential Beltrami operator, as 


Vs 2 (o,,r) — U 52 (z) ■ Ve 3 , 

where V e 3 is the usual Euclidian gradient operator of E 3 and the dot in the r.h.s. denotes 
the 3D tensorial contraction. Note that, obviously, Vs 2 (o,r) = v 52 /R. 

The Laplace-Beltrami operator (or second differential Beltrami operator) is defined in a 
similar way as the restriction of the 3 D Euclidian Laplacian A e 3 to the surface of the sphere 
251. 

A W) - A E3 - - - m . (7) 

We have the scaling relation As 2 ( 0i r) = A s 2 /R 2 and, on the unit sphere, in spherical coor¬ 
dinates 

. 1 d f n d \ Id 2 

A 5 2 = ~ ■ (8) 

sm # 89 \ 89 J sin # dtp 2 

A s 2 identifies with minus the squared angular momentum operator of Quantum Mechanics; 
therefore it is a Hermitian operator the eigenvectors of which are the 3D spherical harmonics 
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y ] m (6,<p) where l is a non-negative integer and the integer m satisfies — l < m < +1. 
Moreover A s 2 ^i n = — l(l + l)y) m . 

A general solution of Laplace-Bcltrami equation in S 2 , i.e. A s 2 Vs 2 = 0, can easily be 
obtained in spherical coordinates by expanding Vs 2 (0,ip) in Fourier series (see e.g. 2c| ) 

m=-\- oo 

v s ,(6, v )= J2 V m (6y mr , ( 9 ) 


which yields for the Fourier coefficients 


sinO-^sinO-^Vye) = m 2 V,a(0). 


( 10 ) 


In the above equation the change of variables x = log tan § (0 < 6 < tt) leads to the 


elementary differential equation 


with the solutions 


= m 2 V m (x ) , 


( 11 ) 


Vo(x) = a + bx , for m = 0 , 
V m (x) = K± m e ±mx , for m/ 0 , 


(12a) 

(12b) 


where a, b , and iF± m are arbitrary complex constants. The general real solution of Laplace- 
Beltrami equation on the sphere can thus be finally written as 


Vs 2 (#> (p) = a 0 + b 0 log tan 


6 tt-A 

n= 1 
oo 

+E 6 

n=l 


tan ■ 


6 ' 

cot- 


e 


cos (rup + a r 
cos (mp + (3 n ) 


(13) 


where a n , b n , a n , and [5 n are arbitrary real constants. The terms of the r.h.s. of (TT3l) 
which are singular at 6 = 0 may be interpreted as the potentials created by point multipoles 
located at the north pole A f and those singular at 9 = tt as the potentials of point multipoles 
located at the south pole S. Remarkably the isotropic term oc log(tan6 l /2) in the r.h.s. is 
singular both at 6 = 0 and 6 = tt and should identify with the potential of a unit bi-charge. 

In order to check this assertion, let us first consider a bi-charge made of a charge +q 
located at point M 0 of ^(O, R) and its companion — q located at the antipodal point M 0 
(with z 0 = —z 0 ). The potential at point M is the solution of Poisson’s equation 


A s 2 V*m 0 ( m ) = -2nq{5 S2 (z 0 , z) - 5 S2 { z 0 , z)} , (14) 
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where the Dirac’s distribution in S- 2 is defined as 5s 2 (z 0 ,z) = 5(1 — z 0 • z) 25]. Eq. (TT4l) is 
solved by expanding both sides on spherical harmonics; after some elementary algebra one 
finds 


07 _i_ I 

V qMo( M ) = £ ' ;( /+ i) ^( COS V’M°m) 

, , fpM 0 M / 1K \ 

= -q - q log tan , (15) 

where Pi(x) is a Legendre polynomial and the prime affixed to the sum in fjl5]l means the 
restriction that l is an odd, positive integer. In this equation we have introduced the geodesic 
distance i/jmqM = arccos(z 0 • z) between the two points M and M 0 on the unit sphere S 2 - 
The first term in the r.h.s. of Eq. (1T3|) is thus indeed the potential created at a source-free 
point M by a point bi-charge located at the north pole Af. 

It is the place to introduce Dirac’s function on a sphere of radius R ^ 1; it will be denoted 
S(Mo,M) = R~ 2 6< s 2 ( z o,z). Since the Laplacian also scales as R~ 2 with the radius of the 
sphere the potential V q ' Mo (M) is independent of the latter. Thus, in (S 2 (0, R), Poisson’s 
equation for a bi-charge takes the form 


A&lo.KVhqV) = —2nq {S(M„, M) - S(M 0 , M)} , 


(16) 


The electric field given by 


ijibi 

^q,Mo 


(M) = --V&1Z“, 0 (M) 
1 1 


-t 


M 0 M 


(M), 


(17) 


R sin V’MoM 

where t m 0 m(M) = cot^MoA / 2 — z 0 /sin denotes the unit vector, tangent to the 

geodesics M 0 M at point M, and orientated from point M 0 towards point M [4,[J Note that 
this vector differs from t m 0 m(M 0 ) = — cot^M 0 Af z o + z/ sin Am 0 m which is the unit vector, 
tangent to the geodesics M 0 M at point M 0 (also orientated from M n to point M). One 


0E|. 


checks that the electric field E b q l Mo {M ) satisfies to Gauss theorem 

Let us consider now a pseudo-charge made of a point charge +q located at at point M 0 
of and a uniform neutralizing background of charge density — g/(47r). The potential 
VqM Q (M) at point M is the solution of Poisson’s equation 

1 


a s 2 V*m o (M) = -27rg{5 52 (z 0 ,z) - —} . 


47T 


(18) 
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Eq. (ITS]) can also be solved by expanding both sides on the complete basis set of spherical 
harmonics with the result 


< ~ X ” > 07 _i_ 1 


2 ; = 1 ^ + 1) 

9 , • i’MoM 

= - q log sm- 

2 y & 2 


(19) 


from which the electric held is readily obtained 


E 


5*0 ( M ) = 


q , i>M 0 M , /, , n 

= 7^ cot —t m 0 m{M) , 


( 20 ) 


AS- 


an expression which can also be obtained by applying Gauss theorem |4, [51. Some remarks 
are in order. 


• (i) One checks that, of course one has : V^ Mq (M) = V^ s Mq (M) — Vki j 0 (M), since the 
backgrounds of the two pseudo-charges canccll each other. 

• (i) For R 1 one has 

A& ( o,B)rr Mo (M) = -2 I r,{i(M„, M) - 1} . (21) 

where S = 47ri? 2 is the surface of ^(0, R ). 

• (iii) Since the 3 D distance between points Mq and M, i.e. the length of the chord 
joining the two points is d = 2Rsm('ipM 0 M/2), it turns out that the expression of the 
potential of a pseudo-charge in <S 2 coincides with that of a point charge in the plane 

e 2 . 


• (iv) Although the potential of a pseudo-charge does not satisfies to Laplace-Beltrami 
equation at source-free points, the potentials created by neutral multipoles made of 
pseudo-charges indeed do, since the backgrounds cancell, as it will be seen below. 


In order to extend the Green’s function expansion (141) to the sphere we have found the 
following trick. Let us introduce, as in Ref. 27], the Cayley-Klcin parameters 


a = exp(i<p/2) cos - , 

Q 

j3 = — i exp(— i(p/2) sin - . 
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(22a) 

(22b) 









A short calculation reveals that, for two unit vectors z\ and Z 2 of d> 2 , one has 


\ct1P2 -CH2AI = sin ^ > 


from which it follows that 


log sin -y- = log sin -y + log cos -y + log 11 — z exp(?A(^) | , 


( 23 ) 


(24) 


where 9 > = sup(0i,0 2 ), 9< = inf(0i,0 2 )j A 93 = ±((pi — ^ 2 ), and z = tan(#</ 2 )/tan( 6 l > / 2 ). 
We note that log |1 — zexp(zA<y 3 )| = (l/2)[log(l — zexp(iA<y 3 )) + log(l — 2 exp(—*Ayj))] and 
that, since \z\ < 1, we can use the expansion of the complex logarithm log(l — Z) within its 
circle of convergence, i.e. log(l — Z) = — Z n /n for \Z\ < 1, where Z = zexp(±?Ay>). 

In doing so, we obtain the following expansion : 


,.■012 . 9 > °< 

— log sm -y- = — log sm —-log cos — 


E 

n=l 


1 


tan(6 l < /2) 

tan(0>/2) 


cos(nAip) , 


(25) 


which seems to be a new mathematical result. In order to make something usefull of that 
esthetic formula, we consider now a set of N pseudo-charges g* of S 2 ( O , R ), with spherical 
coordinates (0j, <£>*), and all around the north pole Af, i.e. ir > 9 0 > 9^ Vi = 1,..., N. A 
multipolar expansion of the electrostatic potential V^ Q R (M) created by these charges at 
some point M = ( 9 , cp) of the surface of the sphere away from A f follows from Eq. (1251) under 
the assumption that n > 9 > 9 0 . One finds 

N N . 

v si(o,R) ( M ) = - J2 f _ J2 ft lQ g sin -y 51 

i=l Z—1 


vr - 5 e <?«- lo s 


Q' 

2 R sin - 
2 


+ EE 


n=1 v 


nX r 


Q n v Gxp(-inutp) , (26) 


where is some unessential constant. In (|26l) the index v assumes, as in Sec. Ill A1 the 
values —1 and +1 only. On the sphere, the complex multipolar moments are defined as 

N 

Qnv = ^2 qiX i ex P {invipi) . (27) 

Z=1 

In Eqs. and (127j) we have introduced the variables X = 2Rtan(9/2) not to be confused 
with the length 2Rsin(9/2) of the chord NM which constitutes the argument of the “log” 
term in the r.h.s. of ([26]) . 
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We note that, since Q n -v = Qnvi there are two independent multipoles of order n only 
(except for the charge of the distribution, i.e. the degenerate case n = 0). In the ther¬ 
modynamic limit p = R6 fixed, R — > oo we have X —> p and the multipole moments (1271) 
as well as the multipolar expansion of the potential in (125]) reduces to their “Euclidian” 
expressions in the plane T{Af) tangent to the sphere at the north pole Af, resp. Eqs. (EJ) 
and ([H]). We remark that quite generally a point multipole located at some point M of S -2 
coincides with a 2 D Euclidian point multipole in the tangent plane T(M) with the same 
complex moments Q nv . Therefore the electrostatic potential Q n; ,exp(— inucp) / (nX n ) in the 
r.li.s. of (1251) may be interpreted as the potential created by a multipole of order n and 
value Q nv located at the north pole Af of the sphere. If n ^ 0 it clearly is a non-isotropic 
solutions of Laplace-Beltrami equation, singular at 9 — 0, cf Eq. ([13]) . 

For the sake of completeness we consider now a set of N bi-charges of S 2 { 0 ,R). Let qt 
(i = 1,..., N) be the point charges located in the northern hemisphere at the points Mj 
with spherical coordinates (#j, </?*). Their N companions —qi are located at the antipodal 
points Mi in the southern-hemisphere, with coordinates &i = n — 6 i and Jp i — ’K + <pi. We 
will denote by V^ 0 R ^ (. M ) the electrostatic potential at some point M of the surface of the 
sphere. Assuming that < 6 < n — 9i, Vi and making use of (125]) one finds 


N 


V£(o,R)( M ) = v S2(o,R)( M ) + + lo S sin 


2=1 


N 


MiM > 


2=1 


v «“ - 


2 R tan 


1 ^ 9 nu ex P (-invip) 


n =1 v 


n[2R] r 


+ 


d 

cot - 

2 


- (-l) r 


e 

tan - 

2 


(28) 


where Eq 111 is some irrelevant constant. The complex multipole moment Q nu has been defined 
in Eq. (1271) . The contribution n = 0 is the potential of a point bi-charge of total charge 
qi and located at the north pole A f. The contribution n ^ 0 is that of a point multipole 
of order n and complex moment Q nv located at Af (and its dumbcll companion at the south 
pole S ). Its electrostatic potential has two singularities in 9 — 0 and 0 = tt (as expected) 
and it is one of the non-isotropic solution of Laplace-Beltrami Eq. (H3l) . To distinguish the 
two types of multipoles encountered in this section, we shall christen mono-multipoles those 
obtained from pseudo-charges and bi-multipoles those obtained from bi-charges. 
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III. DIPOLES 


A. The electric potential of a point dipole 

Henceforth we specialize in the case of dipoles. We start with a point mono-dipole located 
at the north pole A f. It may be seen as a system of N = 2 pseudo-charges of S- 2 ( 0 ,R) : 
a first charge Sq with spherical coordinates ( S9 0 , p 0 ) and a second one with opposite sign 
— 5q and coordinates (69 0 ,po + 7r). We then take the limit 5q —* oo and 59 0 —> 0 with the 
constraint that the dipole moment // = 2 R9 0 5q is fixed. In this limit the vectorial moment 
p = fi (cos <p 0 e x + sin<y 2 0 e y ) belongs to the horizontal plane T(M) = ( e xi e y ) and its two 
non-zero complex multipole moments are Q nu = 5i in /xexp(yz/<p 0 ) where v = ±1. By inserting 
this expression in Eq. (T26|) one obtains the dipolar potential at some point M of the sphere 

V^li°( M ) = 7 ^ cos(cp - ip 0 ) cot Q - . (29) 

In a similar way, one obtains for a bi-dipole located at the north pole M 

v uA M ) = 7^ cos O - Po) | cot ^ + tail 11 • (30) 

Note that Vj^°^°(M) is a solution of Laplace-Beltrami Eq. (TT3T) since the backgrounds of 
the two charges ± 6 q cancell each other. There is a single singularity in 9 = 0. In the 
thermodynamic limit R —> oo, p = R9 fixed one recover the Euclidian dipolar potential 
p/p 2 of a 2D dipole of the (e^, e y ) plane. V$^(M) is also a solution of Laplace-Beltrami 
Eq. (fT5|l but with two singularities at 9 = 0 and 9 = n. Note that the dipoles at the south- 
pole and the north-pole are the same vectors of E 3 . The additional term tan(0/2) in the 
potential vanishes in the thermodynamic limit R —> oo, p = R9 fixed, and does not change 
the conclusion that one also recovers the Euclidian dipolar potential in this limit for the 
bi-dip ole. 

In the general case where the dipole say p 0 of modulus p = ||/^ 0 || is located at some point 
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M 0 of £ 2 ( 0 , R) one obviously have 




/•bi 


2R 

h 


4R sin 2 


/i 1 


R sin M 0 M 

H 1 

R sin 2 V'MoM 


So • t MM 0 {Mo) 


f So • z 

(31a) 

_ ymono / M \ 

M 0 ,p 0 v 

(31b) 

So • 0 ) 


■ So • z , 

(31c) 


where we have introduced the unit vector s 0 = pa/p. 


B. The electric field of a point dipole 


The electric field created by a mono-dipole located at some point M 0 is obtained by 
taking minus the gradient of the potential V^° n ° o (M) at point M with the result : 


KiZS M ) = -R ■ V 5 3 ,m^m 0 °“(M) = 2ttG “(M, Mo) ■ /x 0 , 


where we have introduced, as in Ref. 


(32) 


B. 


the tensorial dipolar Green’s function 


G™™ (M, M 0 ), for which we give two useful expressions : 


G mono 
0 


(M,M 0 


1 


47ri? 2 1 - COS tjj M q m 

[(1 + COSI^Mom) t-MMo(M)t M Mo(M 0 ) — Us 2 (z) • U < 5 3 (z 0 )] , 


1 00 +1 1 


R 2 IQ + 1) 

1=1 m=—l K ’ 


v 53 Fr(z)v 53 b ra ( Z 0), 


(33a) 

(33b) 


as a short algebra will show. We stress that G™ ono (M, M 0 ) is a 3D dyadic tensor of the type 
A(M)A(M 0 ), A(M) and A(M 0 ) being two vectors tangent to the sphere at the points M 
and M 0 , respectively. It is easy to show that in the limit i/>m 0 m —* 0, Go nono (M, M 0 ) tends to 


its Euclidian limit Go,£ 2 (M, M 0 ) = [—' Ue 2 (M) + 2pp\/(2irp 2 ), with p = M 0 M and p = p/p 


and where Ug 2 (M) = e e e d + e^e^, is the unit dyadic tensor in t 


With arguments similar to those exposed in Ref. 17-ll9l. .23] for the 3D case, one can 


re tangent plane 7~(M). 


easily show that the distribution Go ,e 2 (M, M 0 ) has a singularity —(l/2)XJ8 ( ' 2 \p). Therefore 
G™ ono (M, M 0 ) is singular for —> 0, with the same singularity. As for the 3D case [ 7 ] it 
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may be important to extract this singularity and to define a non-singular Green’s function 
GflionoM 0 ) by the relations 


Gg lono (M, M 0 ) = g”’ 5 (M,M 0 ) - -5(M,M 0 )U 52 (z) 


G 


mono, 5 
0 


(M, Mo) 


( G™(M, M 0 ) , for R^m 0 m > # , 

[o , for Rwm 0 m < S , 


(34a) 

(34b) 


where 6 is an arbitrary small cut-off ultimately set to zero. It must be understood that any 
integral involving Q™ 0110 ’ 5 must be calculated with 5 0 0 and then taking the limit 5 —)• 0. 
These matters are discussed in appendix [A] together with some other useful mathematical 
properties of G“ ono (M, M 0 ). 

For the bi-dipoles one finds 


G(f (M, M 0 ) =G“(M, M 0 ) + G“(M, M 0 ) 

= 97r 1 p 2 • 2 1 - (35a) 

2i\R sm 0 Mo m 

[2 cosiIjmoM^mm 0 (M)^mmo(M 0 ) — Us 2 (z) • Us 3 (z 0 )] , (35b) 

= JTj-bjT V S>*7 (z) V 5 ,r, m (zo) • (35c) 

l odd m=—l ' 

in addition to a singularity for M —y Mo (the same as that of G“ ono (M, Mo)) the Green’s 
function Gq‘(M, M 0 ) bears another singularity for M —> Mq. 


C. Dipole-dipole interaction 

We end this section by defining the interaction of two dipoles (M 1; /x 1 ) and (M 2 , /x 2 ) as 
W Mlj/Lt2 = —[i 1 ■ 27tGo(1, 2 ) ■ /x 2 which gives for the mono dipoles with the help of Eq. (13 3a ll 

= ^2 ^ _ cos ^ S 1 ■ S 2 — (1 + 008-012) (tl2(l) • Si)(ti 2 (2) • s 2 )^ , (36a) 

= 7 ^--, -— -r- f s i ' s 2 + t- l —— (si • z 2 )(s 2 • zfi) , (36b) 

2R 1 1 — COS 012 \ 1 — COS 012 ) 

were we have adopted the simplified notation t 12 (0 = (i = 1,2) to denote the 

two unit vectors tangent to the geodesic Mj M 2 . The interaction energy of two bi-dipoles is 
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given by 





(37b) 


(37a) 


One recovers the well-known Euclidian limit r\j ~ (a^ 2 /P 12 )[®i ' s 2 - 2(si • 

p 12 )(s 2 • P 12 )] of the dipole-dipole interaction when ?/q 2 —* 0. 

IV. THE 2D POLAR FLUID IN S 2 (0,R) . 

A. Two equivalent models for the dipolar hard sphere fluid in 52(0,12) 


A fluid of point dipoles is not stable and additional repulsive short range pair-potentials 


must be included in the model to ensure the existence of a thermodynamic limit. In this 
paper we consider only hard core repulsions and thus the DHS model. Two versions are 


possible in S 2 (0,R) depending on whether mono- or bi-dipoles are involved [7]. 

1. Monodipoles 

In this version the mono-dipoles are embedded in the center of hard disks lying on the 
surface of the sphere. I 11 a given configuration, N point-dipoles /x, are located at points 
OM, = Rzi (i — 1,..., N) of 5 2 (0, R ) and the configurational energy reads 



(38) 


where 110 (-0^) is the hard-core pair potential in 5 2 (0,i?) defined by 



(39) 


and is the energy of a pair of mono-dipoles given at Eq. (J36]). Note that the distance 

a is measured along the geodesics and not in the Euclidian space E :i . The hard disks are in 
fact curved objects, similar to contact lenses at the surface of an eye-ball. 


A thermodynamic state of this model is characterized by a density p* = Na 2 /S where 
S = A-kR 2 is the 2 D surface of the sphere S 2 (0,R ) and a reduced dipole moment /T with 
p* 2 = p 2 /(ksTc r 2 ) (ks Boltzmann constant, T absolute temperature). 
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2. Bi-dipoles 


In the second version we consider a fluid of dipolar dumbells confined on the surface of the 
sphere S 2 (0, R ). Both dipoles of the dumbell are embedded at the center of a hard sphere of 
diameter a. In a given configuration, N dipoles /x* are located at points OM; = Rzi and their 
N companions p, at the antipodal points OM; = — (i = 1,..., N) of S 2 (0,R). Each 
pair constitute a bi-dipole and only one of the two dipoles lies in the Nothern hemisphere 
S 2 (0,R) + . The configurational energy reads 

^ N . 1 N 

u ({z is ^}) = v ns(Aj) + ^ ’ ( 4 °) 

*7 

where ^hs (V’ij) is hard-core pair potential defined by 

, 00 if (j/R > iba or iba > tt — a/R , 

4 ( 4 )= ( 41 ) 

I 0 otherwise , 

and the dipole-dipole interaction is that defined at Eq. (]37lh In Eq. (I40li the vectors 

z i can allways be chosen in the northern hemisphere S 2 (0,R) + because of the special sym¬ 
metries of the interaction. Indeed it we have by construction V^. ^ ( M ) = I/A 1 (M) (cf 
Eq. USED) and EA (M) = EA (M) (cf Eq. (!35aD ) for any point M of S 2 (0, R) + and any 
configuration of bi-dipoles. It is thus clear that the genuine domain occupied by the fluid is 
the northern hemisphere S 2 (0,R) + rather than the whole hypersphere. The interpretation 
of the model is therefore the following. When a dipole /x- quits S 2 (0, R) + at some point Mi 
of the equator the same m* reenters at the antipodal point Mi. Therefore bi-dipolcs living 
on the whole sphere are equivalent to mono-dipoles living on the northern hemisphere but 
with special boundary conditions ensuring homogeneity and isotropy at equilibrium (in the 
case of a fluid phase). 

A thermodynamic state of this model is now characterized by a density p* = Na 2 /S 
where S = 2nR 2 is the 2 D surface of the northern hemisphere S^iO^R) and the reduced 
dipole / 1 * with p * 2 = p 2 /(fc^Tcr 2 ) as in Sec. IIV A 11 

B. Thermodynamic and structure 

Most properties of the DHS fluid at thermal equilibrium can be obtained from the knowl¬ 
edge of the one- and two-body correlation functions defined on the sphere and in the canon- 
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ical ensemble as 


p w ( i) 


P (2) ( 1,2) 


N 


Mi) 5 {oli — ati) 


i 2=1 

N N 


^(! - s a) s ( M ^ M i) S(M 2 , Mj) s( ai 

i= 1 3=1 


5(a2 



(42a) 

(42b) 


Here (i) = (z i; a*) denotes both the position z* of dipole “i” and its orientation, specified 
by the angle a*. This angle can be measured for instance in the local basis (eg., e^.) of 
spherical coordinates, i.e. a, = (e^jSj). For a homogeneous and isotropic fluid one has 
p^(l) = p/(27r) where p = N/S is the number density and 

P <2) (1.2)= (^) 2 9(1,2), (43) 


where the correlation function g( 1, 2) has been normalized so as to tend to 1 at large dis¬ 
tances. To exploit the symmetries of the fluid phase It is far more convenient to use, instead 
of the ex*, the angles fa = (ti 2 (i),Sj) (i = 1,2). In order to complete the local basis in the 
tangent planes % (i — 1, 2) one defines ni 2 (i) = z* x ti 2 (i) (i = 1, 2). Remarkably the vector 
ni 2 is invariant along the geodesic MiM 2 and one has 


ni 2 m = n 12 = 


Zi X z 2 


(44) 


sin -0 1 2 

Note by passing that the unit dyadic tensor U < s 2 (z i ) = ni 2 ni 2 + ti 2 (i)ti 2 (i) for i — 1,2. 

These variables (fa, fa) serve mainly to introduce a basis of 2D rotational invariants 
'f?rnn(/3 1 , fa) upon which to expand the pair correlation function g( 1,2). The rotational 
invariants were introduced by Blum and Toruella in the context of 3 D homogeneous and 


isotropic molecular fluids 


and their 2D counterparts in the 2D plane were given in 


Refs [4, 5]. On the sphere the extension of the latter result is straightforward and these 
invariants read 


^«m(l> 2 ) = ex P (imfii + m/3 2 ) . 


(45) 


These functions are indeed invariant under the action of a global rotation of E 3 of center O. 
Moreover these invariants are orthogonal is the sense that 

< Vmn^pq >= J ^ J ^T mn (l, 2)^(1, 2) = 5 mp 5 nq . (46) 

and in an homogeneous and isotropic phase of the DHS fluid one has the expansion 

OO 

5 , ( 1 ’ 2 ) = 5Z 9™n(r 12 )^ mn(l fa) , (47) 

m,n =—00 
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where the coefficients g mn (j' 12 ) are functions of the length r 12 = Ri/j 12 of the geodesic MiM 2 . 

Reality of the correlation function g( 1, 2) and invariance through the reflexion ((3i, (3 2 ) -H- 
(— Pi, — /3 2 ) impose that g mn = g- m - n = Qmn- These conditions enable us to define the real 
rotational invariants 

&mn = g (Tmn “I - ^—m—ri) 

= cos (m/3 1 + n/3 2 ) . (48) 


To give some physical and geometrical interpretation of some of these invariants we first 
recall the expression of the scalar product on the sphere 1 S n (0,R) [4-7]. Taking the scalar 
product A(l,2) of two vectors Si and s 2 located at two distinct points, Mi and M 2 of 
1 S 2 (0,i?) needs some caution. It first requires to perform a parallel transport of vector S[ 
from Mi to M 2 along the geodesic MiM 2 and then to take a 2D scalar product in space 
T(M 2 ). Thus 

A(l, 2) = Ti 2 S! • s 2 , (49) 


where, in the r.h.s. the dot denotes the usual scalar product of the Euclidian space T(M 2 ). 
Vector ri 2 Si results from a transport of Si from the space T(Mi) to the space T(M 2 ) along 
the geodesic MiM 2 , keeping its angle with the tangent to the geodesic constant. Explicitely 
one has: 

T 12 S 1 = Si - Sl 22 (zi + z 2 ) , (50) 

1 + cos i>i 2 

and thus 

A(l,2) = Sl -s 2 - (Sl ^ 2)(S2 ' Zl) . (51) 

The above expression (j5Tj) is in fact valid for a hypersphere S n (0, R) of arbitrary dimen¬ 
sions. In the case n — 2 more simple expressions can easily be obtained with the help of 
the angles /% introduced previously. Starting from Si = cos / 5 1 t 12 (l) + sin / d 1 n 12 one deduces 
from ([50]) that r 12 Si = cos / 0 1 t 12 (2) + sin/3in 12 from which it follows that 


A(l,2) = cos(/?i — fa) , 


(52) 


and therefore A = 4 ) i_i. 

By analogy with the 2D Euclidian case we also introduce the manifestly rotationally 
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invariant function .0(1,2) 


0(1, 2) = = 2( Sl • t 12 (l)) (s 2 • t 12 (2)) - A(l, 2) 
_ _ (si • z 2 )(s 2 • Zi) 

1 2 1-COS-012 

An explicit calculation shows that in «S 2 (0, R) one also has 

0(1,2) = cos(£i + p 2 ) , 


(53) 


(54) 


and therefore O = $n. 

The angular Dependance of the dipole-dipole interactions (136|) and (137)) should be rota- 
tionally invariant and indeed one finds that it is a combination of the invariants O and A. 
More precisely one has 


Titmono 


b 


1 + COS 012 


w, 


bi 

Ul,M2 


R 2 sin 2 0i2 
b 2 

R 2 sin 2 0i2 


0 ( 1 , 2 ) 


1 + cos 


012 A 




0 ( 1 , 2 )- 


1 — cos 


2 —) A(l,2) 


(55a) 

(55b) 


The expressions of the projection g mn (r i2 ) of the pair correlation g( 1, 2) on the rotational 
invariants are easily deduced from the definition (I42bj) of the two point function p^(l,2) 
and the properties of orthogonality of the $ mn . One finds (in the canonical ensemble) 

1 1 / ^rnn (0j)x(0ij -0 12 ) 


g mn (r 12) = 




mnl^mn) Np \ ‘ZliR sill(0ij)50 


(56) 


where 50 is the bin size and x i s defined as 


X '(0 - 012 ) 


1 if 012 < 0 < 012 + 50 
0 otherwise. 


(57) 


Note that ( D\D ) = (A|A) = 1/2 while (<h 00 |cl> 00 ) = 1. We want to point out that for the 
DHS fluid of mono-dipoles 0 < 0i 2 = r^/R < 7 r, however, for the fluid of bi-dipoles, only 
the range 0 < 0 i2 < 7t/2 is available, because of the special boundary conditions involved 
in the model. 

If follows from the precedent developments that, at equilibrium, the mean energy per par 
particles reads 




< y , w mono > 


V 

’4 


2 N 

T 

d0 sin0 


/l D (i?0) 

1 — COS 0 ’ 


(58a) 
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for mono-dipoles and 


/3« bi = 


< £*,- WS > 


4 j 
and 


2iV 

t/2 


dsini/j 


h D (Ri/}) h /1 (R'ijj) 


1 — COS -0 1 + cos-0. 


(59a) 


we have introduced the dimensionless parameter y = 


for bi-dipoles. In Eqs 
7 Tpfdfj? and h = g — l. 

For the sake of exhaustivity we finally quote the expression of the compressibility factor 

~R 


^mono (bi) _ ^ _|_ ^mono (bi) q. dp* 


a 

- Sill — 
<7 it 


<? U > + 0) 


(60) 


Note the prefactor of the isotropic projection g 00 (<J + 0) at contact which accounts for 
curvature effects. The usual Euclidian expression = l + /3uoo + (np*/2)g oo (a+0) emerges 
in the limit R —» oo, both for mono- and bi-dipoles, as it was expected. 


C. Fulton’s theory 

Let us consider quite generally a polar fluid occupying a 2D surface A with boundaries 
<9A. We assume the system at thermal equilibrium in a homogeneous and isotropic fluid 
phase. The fluid behaves macroscopically as a dielectric medium characterized by a dielectric 
constant e. Due to the lack of screening in such fluids, the asymptotic behavior of the pair 
correlation function is long ranged and depends on the geometry of the system, i.e. its 
shape, size, and the properties imposed to the electric held (or potential) on the boundaries 
<9A as well. As a consequence, the expression of the dielectric constant e in terms of the 


fluctuations of polarization also depends on the geometry. T 


taken into account in the framework of Fulton’s theory 171-ll9l| which achieves an elegant 


rese issues can be formally 


synthesis between the linear response theory and the electrodynamics of continuous media. 

This formalism can be extended without more ado to non-euclidian geometries and was 
applied notably for 3D hyperspheres and cubico-periodical systems in Refs. I?, Qj- In this 
section it is applied to the 2D sphere both for mono- and bi-dipoles (a separate treatment 
of each model is however necessary). To our knowledge the case of the 2D plane E 2 was 
never considered in the framework of Fulton’s theory and it is given a special treatment in 
appendix 
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1. Mono-dipoles 


We consider a fluid of N mono-dipoles in S2(0,R ) at thermal equilibrium in the pres¬ 
ence of an external electrostatic field £{M) E T(M). The medium acquires a macroscopic 
polarization 

P(M) =< P(M) > f , (61) 

where the brackets denote the equilibrium average of the microscopic polarization P(M) 
in the presence_of the external field £. In tS 2 (0,i?) the microscopic polarization P(M) is 


defined as 


7,l2S| 


N 


P(M) = ^U 52 (z)-^. 5(M,M j ), 

3 = 1 


(62) 

where the unit dyadic tensor Us 2 (z) in the r.li.s. of Eq. (j521) ensures that vector P(M) 
belongs to the plane T(M) tangent to the sphere at point M. 

The relation between the macroscopic polarization P(M) and the external field £(r) 
can be established in the framework of linear-response theory, provided that £ (M) is small 
enough, with the result 


2nP(M) = X° £ (= [ 
V Js 2 


dS' x(M,M')-£(M' 


(63) 


>S 2 (0,R) 

The r.h.s. of Eq. (T63|) has been formulated in a compact, albeit convenient notation that 
will be adopted henceforth, where the symbol o means both a tensorial contraction (denoted 
by the dot ” • ”) and a spacial convolution over the whole sphere. 

From standard linear response theory the tensorial susceptibility x is then given by 


X{M U M 2 ) = 2tt/3 < P(M!)P(M 2 ) > , 


(64) 


where the thermal average < ... > in the r.h.s. of (RhO) is now evaluated in the absence of 
the external field, £ = 0. The susceptibility tensor x(M 1? M 2 ) can be expressed in terms of 
the one- and two- point correlation functions (1421) as 


%(Mi, M 2 ) — x 5 (Mi, M 2 ) + x d (M : i, M 2 ) , 


(65) 


where the “self” part reads 


Xs(M 1 ,M 2 )=y6{M 1 ,M 2 ) U*( Zl ) , 


( 66 ) 
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and 

xd(m 1 ,m 2 ) = 2 (67) 
After making use of the expansion (1471) of h( 1, 2) in rotational invariants one finds finally 

X (M 1 ,M 2 ) = yS(M 1 ,M 2 ) U* ( Zl ) 

+ Y { [(! - cos "012)^ A (?" 12) + (1 + cos^i 2 )h D (ri 2 )] ti 2 (l)ti 2 (2) 

+ (h A (r 12 ) - h D (r u )) U l 5 2 (z 1 ) • U 52 (z 2 )} . (68 ) 

We point out that in the limit r i2 fixed, R —> oo the above expressions reduces to that 
derived in Appendix iBl for the 2D euclidian space (cf Eq. (1B6|) ). 

The dielectric properties of the fluid are characterized by the dielectric tensor e which 
enters the constitutive relation 

27tP = (e — I) o E , (69) 

where E denotes the Maxwell field and I(Mi,M 2 ) = Us 2 (zi)5(Mi, M 2 ). The Maxwell field 
E(M) is the sum of the external field £(M) and the electric field created by the macroscopic 
polarization of the fluid. Therefore one has 


E = £ + 2ttG“ o P . (70) 

It is generally assumed that e is a local function, i.e. e = el. More precisely, it is plausible 
-and we shall take it for granted- that e(Mi,M 2 ) is a short range function of the distance 
between the two points (Mi,M 2 ), at least for a homogeneous liquid, and one then defines 

eU 52 ( Zl )= f dS 2 e{M 1 M 2 ). (71) 

Js 2 {o,r) 

In general (e — I) ^ x since the Maxwell field E(M) and the external field £{M) do not 



E = £ + G mono o <r o £ , (72) 


where cr = e — I and G mono (Mi, M 2 ) is the macroscopic dielectric Green’s function defined 
by the identity 


G mono mono ,'i _ /—(mono 

= G 0 O (1 — <7 O (jr Q 


,-l 


(73) 


where the inverse must be understood in the sense of operators. It is easy to show that the 
electric field created at a point M x e S 2 (0, R ) by a point mono-dipole /x 2 located at M 2 in 
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the presence of the dielectric medium is then given by 27rG mono (Mi, M 2 ) ■ /i 2 . This remark 
enlightens the physical meaning of the macroscopic Green’s function. 

Combining Eqs. (1551) . (1551) . and (1751) yields Fulton’s relation 

X = cr + cr o G mono o cr . (74) 

To go further one has to compute the macroscopic Green’s function G mono . Our starting 
point is the following identity, proved in appendix [7] : 


G mono pmono _ mono 

0 0 ^0 — -,or 0 


(75) 


Therefore —G™ ono is a projector and has no inverse. Assuming the locality of cr one is then 
led to search the inverse (I — cr o G™ ono ) -1 in the r.h.s. of (175]) under the form al + bG. 


mono 

0 


where a and b are numbers (or local operators). By identification one finds a = 1 and 
b = er/(l + a) yielding for G mono the simple (and expected) expression 


G 1 


G mono //-i 1 _\ _ nmono /_ 

0 /( 1 + M = G 0 /e. 


(76) 


The results derived above allow us to recast Fulton’s relation (175]) under its final form 

^l ) 2 


1)1 (M 1 ,M 2 ) + 


Gr n °(Mi,M 2 ) = X (M 1 ,M 2 ) 


(77) 


We stress that the above equation has been obtained under the assumption of the locality 
of the dielectric tensor e(Mi, M 2 ). Therefore it should be valid only asymptotically, i.e. for 
points (Mi, M 2 ) at a mutual distance ri 2 larger then the range £ of e(Mi, M 2 ). 


Expressions for 


Following Refs. 


B 


he dielectric constant are obtained from Eq. (177]) by space integration. 


2a 


30] one integrates both sides of Eq. (177]) and then takes the trace. 


The integration of M 2 is performed over a cone of axis z 1 and aperture £; 0 and then M x is 
integrated over the whole sphere S 2 (0,R). The singularity of the dipolar Green’s function 
G 0 (Mi,M 2 ) for f; 12 0 must be carefully taken into account and this delicate point is 

detailed in appendix [A] (see Eq. (1A9j) ). One hnds hnally 


i+k 


i)= 


4e 


(1 + COsf> 0 ) = m2 (V’o) With 0 < ifo < 7T , 


(78) 


where the dipolar fluctuation m 2 (£7 0 ) reads as 


m 2 ( , 0o) = 


7 T/3fi 
~S 


2 N N 

0 (Vh - Aj) > , 

i j 


(79) 
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where Q(x) is the Heaviside step-function (Q(x) = 0 for x < 0 and Q(x) — 1 for x > 0). 

We have thus obtained a family of formula depending on parameter 0 < i/jq < 7r; clearly 
they should be valid only if Rij ) 0 is large when compared to the range £ of the dielectric 
constant. The numerical results of Sec. (jVJ) show that this range is of the order of a few 
atomic diameters. It is also important to note that for = tt Eq. (178|) involves the 
fluctuations of the total 3D dipole moment of the system. However, the resulting formula 
i.e. (e — l)/e — itl 2 (7t) = tt(3 < M 2 > /S (with M = pi, the total 3D dipole moment of 
the sphere), albeit simple, is not adapted for numerical applications since, for large values 
of the dielectric constant, a reasonable numerical error on e requires a determination of 
m 2 (7r) with an impractical precision. The choice i(jq = 7t/2 yields the less simple formula 
(e — l)(e + 3)/e = 4m 2 (7r/2) which however allows, by contrast, a precise determination of 
e. Indeed, let 6 e be the error on e, then, for high values of the dielectric constant the errors 
on m 2 (7r/2) and e are roughly proportional, i.e. 5e ~ 4 4m 2 (7r/2). 

The fluctuation m 2 (^ 0 ) is related to the so-called Kirkwood factor G K (i/jq). according to 
the relation m 2 (^o) = yG K ('ip l o) with 

ri’ o 

G K (r ) = 1 + pvr R 2 d^ sin £; h Kirk (i?^) , (80) 

Jo 

where 

^Kirk.^ _ 1 (i cos -0)h A (r) — i (1 — cos ijj)h D (r) . (81) 

Fulton’s relation (177]) can also be used to determine the asymptotic behavior of the 
projections h A (r ) and h D {r ) the pair-correlation function g{ 1,2). By comparing Eqs. (I08D 
and (1771) one obtains readily that, asymptotically, i.e. for large r = Ri/j and ^< 77 , one has 


h 

h 

h 


A 

asymp 

D 

asymp 

Kirk. 

asymp 


Jr) ~ 0 , 


( r ) ~ 
J r ) ~ 


( e — i ) 2 i_L_ 

ype 2 itR 2 1 — cos 

(^-l) 2 1 

ype 47t R 2 


^ ’ 


(82a) 

(82b) 

(82c) 


We note that h A (r ) has no tail and is thus a short-range function. It thus presents the same 
behavior as in the 2D Euclidean space (see Refs. 4, |5] and the new analysis in appendix |B] 
also based on Fulton’s relation). We also note that, in the thermodynamic limit R —> oo and 
with r>( fixed but large, one recovers the expected Euclidian behavior h£ ymp (r) ~ (e — 
l) 2 /(Trype) x 1/r 2 valid for an infinite system without boundaries at inhnity( see e.g. [4s, 5] and 
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Eq. (1B12D of appendix [B]). The Kirkwood’s function /i Kirk ' (r) exhibits a constant tail which 
tends to zero as the inverse of the surface S = AnR 2 of the system. However the integration 
of this tail over the volume gives a non-zero contribution to the dielectric constant. In the 
limit r = R^> £ one has 


G k Oo) = G“ + vrpi? 2 I dij) sin ^ (ity) 


iK 


cipo 


Kirk. 


= Gx- 


1) 


asymp. 


4 ye 


where G 1 ' is the Euclidian Kirkwood factor 


Gl = 1 + 


P 


(1 - COS^o) , 


dr2nrh^(r) , 


(83) 


(84) 


where we have noted that, in the thermodynamic limit, r fixed and R —> oo, /i Kirk '(r) = 
h^(r). Reporting this asymptotic expression of the Kirkwood factor in Eq. d79|) one recovers 
the formula valid for the Euclidian plane which reads as 

( C -l)(e + l) 


2e 


= yG 


K 

OO 


(85) 


which coincides with the Eq. (IB 10D obtained in the appendix [B] for the infinite plane E 2 ] 
note that Eq. [85] can be derived by a whole host of other methods, cf Refs BQ. 


2. Bi-dipoles 

The passage from mono- to bi-dipoles requires some adjustments. First, the appropriate 
bare Green’s function is now Gg 1 as defined at Eqs. f[35]h Secondly the operation o defined 
at Eq. (1631) now involves a spacial integration the support of which is only the northern 
hemisphere 5^(0, R) rather than the entire surface of the sphere. For instance, as shown 
in appendix lAl one has 

[G^ O Gg 1 ] (Mi, M 2 ) = / dS Go i (M 1 , M) • G h 0 \M, M 2 ) 

JS+(0,R) 

= -G^ i (M 1 ,M 2 ) (86) 

for points points and M 2 belonging to the northern hemisphere. Therefore, reasoning as 
in Sec. II V C 11 we find for the dressed Green’s function, asymptotically 

G bi = G bi o (I - <r o G bi ) _1 = Gg’/e . (87) 
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Clearly the susceptibility tensor % D (Afj, M 2 ) keeps the same expressions (1641) and (1681) with 
the restriction that the two points Mi and M 2 both belong to the northern hemisphere 
S^iP^R). As in Sec. IIV C II the dielectric constant can be obtained from Fulton’s relation 
a o G bl o CT and reads now 


A = cr 


e — 1 


(e-l) S 

2e 


cos0 o = m ' 2 fV^o) with 0 < 0 O < 7r/2 , 


( 88 ) 


where the fluctuation m 2 (0 o ) is still given by 


m 2 ('0o) = 


7T/3yU 2 

~s~ 


N N 


<^2^2 S i- S i 0(^0 - i>ij) >, 


(89) 


but with S = 2nR 2 (surface of the northern hemisphere) and 0 < 0 O < 7 t/ 2. The delicate 
mathematical integration of the Green’s function G bl required to derive Eq. (1891) is explained 
in appendix lAl 

We can note that for 0 O — 7r/2, i.e. when the fluctuation of the total dielectric moment 
of the sample is taken into account one still has (e — l)/e = ir/3 < M 2 > /S with M = 
the total 3 D dipole moment of the northern hemisphere. Quite remarkably one obtains 
for 0o = 7r/3 the relation (e — l)(e + 3)/e = 4m 2 (7r/3) which is formally identical to that 
obtained for mono-dipoles with the choice 0 O = 7r/2. It is not a surprise since, in both cases, 
m 2 (0o) accounts for the dipole fluctuation of half the total domain available to the dipoles 
of the system. 

The asymptotic behavior of the pair correlation function is obtained in the same vein as 
in Sec. IIV C 11 i.e. for large r = Ri(j and ^ < 7r/2, one has 

(e-1) 2 1 1 


h 


asymp 


, D 


.(0 


h 1 ' ( r) ~ 

' ^asymp. V' ) 


ype 2 ttR 2 1 + cos -0 ’ 

-l) 2 1 1 

COS 0 ’ 


h Kirk 


asymp 


P) 


ype 2 ttR 2 1 
(e-1) 2 1 


(90a) 

(90b) 

(90c) 


ype 2nR 2 

In the limit R —> 00 and with r £ fixed but large, one recovers once again the 
expected Euclidian behavior h£ ymp (r) ~ (e — l) 2 /(7r ype) x 1/r 2 . By contrast, in the same 
limit, one obtains that /i0 symp (r) ~ —(e — l) 2 /(47 Type) x 1/R 2 which tends to zero for the 
infinite system for which R —> 00 . This behavior is in agreement with the expected short 
range behavior of the projection h A (r) in the 2D infinite Euclidian plane (cf Refs. 0,0 


26 















and the appendix [B]). The Kirkwood’s function /i^* p (r) presents the same constant tail 
— (e — 1 ) 2 /(ype) x 1/S as in the mono-dipole cases. Finally, in the thermodynamic limit, 
the Kirkwood’s factor behaves as 0 ) = m2 (^o )/y — — (e — l) 2 (1 — cos'i/> 0 )/(2ye ) ; 

by inserting this expression in Eq. (jggj) one recovers, as for mono-dipoles, the formula (j85|) 
valid for the Euclidian plane E 2 . 


V. SIMULATIONS 


The only published MC simulations of the 2D DHS fluid the author is aware of are 
those of Morriss and Perrarn 21] published 30 years ago. Their simulation cell is a square 


with periodic boundary conditions and the authors make a correct use of Ewald dipolar 
potentials (an alternative version of that used in their paper, however non implemented in 
actual numerical simulations, is discussed in Ref. 22]). We have retained one of their points, 
i.e. our state I : (p* — 0.7, /i* 2 = 2) as a benchmark. We also report MC data for a second 
reference state, our state II = (p* = 0.6, p* 2 = 4) characterized by a much higher dielectric 
constant. Both states undoubtly belongs to a stable fluid phase as it will be discussed below. 

We performed standard MC simulations of a DHS fluid in the canonical ensemble with 
single particle displacement moves (translation and rotation). Each new configuration is 
thus generated by the trial displacement of a randomly chosen single dipole by means of a 
new algorithm explained in appendix O 

We considered both mono- and bi-dipoles on the sphere and studied finite size effects on 
the thermodynamic, structural and dielectric properties of the fluid. We studied systems 
of N = 250, 500,1000, 2000, and 4000 particles and typically we generated A^conf = 16 10 9 
configurations for each considered state. 

Our results are resumed in Table Q] (state I) and Table ITT! (state II). The errors reported 
in these tables correspond to two standard deviations in a standard statistical block analysis 
of the MC data 31]. 


As apparent in Fig. [T] the reduced internal /3u energies converge linearly with 1/N to 
their thermodynamic limit (duoo for systems involving more than N = 250 particles. This 
is true for systems of mono- and bi-dipoles. This allows a determination of with a 
precision of ~ 3 10~ 5 . The asymptotic values obtained for mono- and bi-dipolcs are identical 
within the error bars. Clearly the convergence towards the thermodynamic limit is faster for 
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bi-dipoles. Probably because, for a given density and number of particles, the radius R of 
the sphere is y/2 larger for a system of bi-dipoles and therefore curvature effects are smaller. 
We note that our finding /3uoo = —1.79002 ± 510~ 5 compares reasonably well with that 
obtained by Morriss and Perram for samples of IV = 100 particles which lie in the interval 
—1.795 < ... < —1.780 depending on the type of boundaries surrounding the system at 
infinity 21]. 

Size effects on the short range part of the projections g 00 (r ), h D (r), and h A (r) are very 
small and cannot be appreciated on the graphs of the functions. To give an idea of these 
effects we give the contact value of these three projections in the tables. These values, 
which result from a Lagrange’s polynomial interpolation, are much less precise than the 
values obtained for the energy. The compressibility factors Z were deduced from Eq. (16(71) 
and are also reported in the tables. The thermodynamic limit of these quantities do not seem 
to obey a simple linear regression with 1/N. However, the low precision on the contact value 
g 00 (u) precludes an unambiguous study of the thermodynamic limit in this case. We can 
infer from our MC data that, for state I one has in the thermodynamic limit = 3.98±0.01 
a value once again not far from that of Ref. 21] Z ~ 4.06 obtained for a small system of a 
N = 100 dipoles. 

The dielectric constants were obtained from Eq. (1751) with t/’o = 7 t/ 2 for mono-dipoles and 
from Eq. (159]) with -0 O = 7t/3 for dipoles. As apparent on Figs. [2] the numerical uncertainties 
on e are large and seem even to increase with system size but do not mask significant finite 
size effects. As for the pressure, it was impossible to deduce from our MC data a convergence 
law of e(N) for N —)■ oo despite the huge number of MC steps involved. However we can 
claim that, for state I one has e^ = 18.0 ± 0.2, which differs significantly from the value 
reported by Morriss and Perram e = 16.0 ± 0.5. 

In order to test the validity of Fulton’s theory one can examine Fig. [3] which displays the 
Kirkwood factor G K (r) as a function of the distance r = R'lp o which defines the radius of 
the cap where the partial electric moment fluctuations are taken into account. As discussed 
in Sec. IIV Cl for r > £ (£ range of the dielectric tensor) this function should be given by the 
asymptotic behaviors given by Eq. (1751) for mono-dipoles and Eq. (1551) for bi-dipoles. The 
dielectric constant used for these asymptotic forms are those given in Tables [T] and [TT1 The 
excellent agreement beetween the MC data for G K {r) and the theoretical predictions even 
yields an estimation for the parameter £, i.e. £ ~ 7cr for both states 1 and II.. 






We check in Fig. [4]that the asymptotic behavior of the projection h A {r) is indeed governed 
by the law of macroscopic electrostatics in <S 2 , he by Eqs. (1821) for mono-dipoles and Eqs. (j90j) 
for bi-dipoles. In all cases for r > £ the asymptotic behavior of h A (r ) is in excellent 
agreement with the theoretical prediction. Note that for mono-dipoles h A (r) has no tail 
as expected from our theoretical developments. In the same vein Fig. 0 displays the ratio 
h D (r) / h® ymp (r) for states I and II and several system sizes. Once again a satisfactory 
agreement simulation/theory can be observed. Note that for mono-dipoles (top curves) the 
noisy behavior at large distances is a consequence of the divergence (due to an insufficient 
sampling) of the term ~ 1/sin -0^ for 7 /;^ ~ n in the l.li.s. of Eq (IMj) which defines h D {r = 
Rif>) as a statistical average. 

The fact that, for both states I and II, the theoretical asymptotic behavior of h A {r) and 
h D (r ) is reproduced by the MC data is a clue that the system is in an isotropic liquid phase. 
At lower temperatures chains of dipoles, rings, and topologically complex arrangements arise 
and the dielectric tensor either does not exist or is no more local. This point will be discussed 
elsewhere. 

VI. CONCLUSION 

In this work we have clarified the laws of electrostatics on the sphere <S 2 (0, R) by introduc¬ 
ing two types of basic elements : pseudo- and bi-charges from which two types of multipoles 
can be obtained : mono- and bi-multipoles respectively. The electrostatic potentials of these 
multipoles can be obtained explicitly and exhaust all possible solutions of Laplace-Beltrami 
equation. 

The application to polar fluids has been discussed in depth. Such a fluid can be repre¬ 
sented either as an assembly of mono-dipoles living in <S 2 (0, R) or as a collection of bi-dipoles 
living in the northern hemisphere iS 2 (0,i?) + . 

Since the dipolar Green’s function are explicitely known for both mono- and bi-dipoles, 
Fulton’s theory of dielectric fluids can be explicitely worked out in the two cases. This 
includes the theory of the dielectric constant of an homogeneous and isotropic fluid and the 
derivation of the asymptotic tails of the projections h A (r) and h D (r ) of the equilibrium pair 
correlation function. It is likely that a violation of these asymptotic laws signals a phase 
transition towards a non-fluid phase. 
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We have reported and discussed MC simulations of two states of the DHS fluid which 
are both in the fluid phase and can serve as benchmark for future numerical works. Fi¬ 
nite size effects have been studied and allow to reach, at least for the internal energy, the 
thermodynamic limit with a high precision (~ 10 -5 ). Such a precision is smaller than the 
precision which could reasonably be obtained on the dipolar Ewald potential and a fortiori 
on the mean internal energy which could be computed by an actual MC simulation which 
would make use of it. This suggests that ^(0, R) is a good geometry for precise MC sim¬ 
ulations of the fluid phase. Simulations involving bi-dipoles appear to attain more quickly 
the thermodynamic limit than simulations involving mono-dipoles. An exploration of the 
phase diagram of the 2D DHS fluid will be presented elsewhere. 
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Appendix A: Mathematical properties of the dipolar Green’s functions in S 2 

1. Convolutions 


First we consider the case of mono-dipoles. Let zi and Z 2 be two points of the sphere S 2 , 
we will show that 


G mono _ /'-'imono/_ _ \ / jn/r, \ mono/_ _ \ pmono/_ _ \ 

0 0 v *0 ( z l 7 z 2 )= / «h(z 3 ) Gq (Zi,Z 3 )-G 0 (z 3 ,z 2 ) 


's 2 


G monop _ \ 

0 V Z 1 1 z 2 ) j 


(Al) 


which shows that —G™ ono is a projector and has thus no inverse with respect to the binary 
operator “ o ”. In order to prove Eq. (lAlj) we rewrite the expansion (133bD of the dipolar 
Green’s function in spherical harmonics as 

OO 

Gr n >i.z 2 )=5]G'(z 1 ,z 2 ), (A2) 

1=1 


with 


Gofa.zj) 


E V & F7(zi)V & Vr(z2) . 

' ' m=—l 


(A3) 
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To make further progress we need apply Green-Beltrami theorem in 5 2 which states 


that 


25j : 


[ dn\7 s J-X7 S2 g = - [ dDfA S3 g, (A4) 

J S2 J S2 

where /(z) and g( z) are functions defined on the unit sphere <S 2 . The proof of theorem (jA 
is not so difficult and can be found, e.g. in the recent textbook by Atkinson and Han 
It follows from Eq. (1A4D and the properties of spherical harmonics that 


25] 


]Gq o Gq](zi, z 2 ) — —Si }P Gq(zi, z 2 ) , 


(A5) 


from which the identity (lAlj) is readily obtained. 

The case of bi-dipoles is slightly different and has already been discussed in Ref. [7] for 
the 3D hyper-sphere £ 3 . We now consider two points Zi and z 2 of the northern hemisphere 
5^. We will show that, again 

Gq 1 o Go‘(zi, z 2 ) = [ dfl(z 3 ) Go^znzs) ■ Gq 1 (z 3 ,z 2 ) 

Js+ 

= “G{j 1 (z 1 ,z 2 ) . (A 6 ) 


Note that, in this case, the space integration is restricted to the northern hemisphere. We 
follow the same strategy as for mono-dipoles and start with the expansion (I35cji of the 
dipolar Green’s function Gq 1 (z 1 ,z 2 ) in spherical harmonics 

G^(zi,z 2 ) = 2^G[ ) (z 1 ,z 2 ) . (A7) 

l odd 

For an odd value of l we have = —l^ m (z) which implies that 

Gj, o GJ(z„ z 2 ) = -i 6,0 G‘ (Zi,z 2 ) , (A8) 

from which Eq. (1A6D follows. 

In this Sec. IA II we implicitely assumed that R — 1 . The reassessment of Eqs. (IA1|) 
and (IA 6 D in the case R 7 ^ 1 is however trivial since the mono and bi dipolar Green’s function 
G 0 (1,2) scales as R~ 2 with the radius of the sphere. Clearly Eqs. (IA1|) and (IA 6 D remain 
valid for R 7 ^ 1 with the replacement dfl(z) —>- dS(M). 
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2. Integration over cones 


This section is devoted the integration of G™ ono (zi, Z 2 ) on a cone of axis zi and aperture 
0 < ipo < We shall prove that 


A CM ^ \ pmono^ „ _ COS 3 

d\l{ Z 2 ) G 0 (Zi,Z 2 ) — - --U5 2 (ZiJ 


(A9) 


7o<V'i2<ho 

where ^12 = cos _1 (zi ■ z 2 ) is the angle between the two unit vectors z\ and z 2 . To prove 
Eq. (lA9j) one needs to take some precaution because of the singularity of G“ ono (z 1 , z 2 ) at 
ip \2 —> 0. We make use of the decomposition (I34al) to rewrite 

[ dSl( z 2 ) G“(zi,z 2 ) = ~U 5a ( Zl ) + lim [ <m( z 2 ) G“(zi,z 2 ) • (A10) 

J 0<pl2<V>0 Z 5 ^° JS<1pl2<i >0 

The integral tf° in the r.h.s. of Eq. (IA10D is computed by using spherical coordinates to 
reexpress the formula (133bjl of the Green function and performing explicitely the integrals. 
One has 


0 (Zl,Z 2 ) = — 


1 f 1 + COS 


ti 2 (zi)t 12 (z 2 ) 


U 52 ( Zl )-U 52 (z 2 


47T ( 1 — COS Ip 12 1 — COS "012 

By taking the north pole at point zi we have the identifications -0i 2 —* 9, ti 2 (z 2 ) —>■ eg(9, p) 
and t 12 (z!) —> e e (9 = 0, ip) = (cos p, sin <p, 0) 1 . Therefore 


jho _ 1 

15 “47T 


/»27T 

sin QdQ / dp 


1 + cos 6 
1 — cos 6 


e e {6,p)e e (Q,p) 


1 — cos 6 
cos Vh — cos 5 


U 52 (0,^)-U 3,(9, p) 
U 52 (zi) • 


(All) 


Taking the limit S —> 0 and gathering the intermediate results (1A10D and (1A11I) one obtains 
the desired result (lA9f) . 

The same type of calculation yields, for bi-dipoles, 


dtt(z 2 ) GS'(zi,z 2 ) = - l)U & (z,) 


(A12) 


J 0<ll'12<1pO 

with the restriction that 0 < ipo < n/2. 
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Appendix B: Fulton’s theory in the plane. 


In this section the domain occupied by the fluid consists of the entire 2D Euclidian plane 
E 2 with no boundaries at infinity. We consider Fulton relation (Thill in this geometry. The 
susceptibility tensor x(Mi,M 2 ) at thermal equilibrium is 


x(Mi, M 2 ) = 2 t t/3/i 2 (^2 s i s j ( r i - r *) ^ (2) ( r 2 - r j)) , 

= %5( r 12) +Xd(*12) , 


where 


OaI; = r, = 


J^ii Vi 


and r 12 = r 2 — iq. The self-part Xs i n HBlbjl reads 


(Bla) 

(Bib) 


Xs(ri 2 ) = Z/I(ri,r 2 ) , 


(B2) 


with y = it p/3 p 2 and I(r l5 r 2 ) = Uh^ 2) (r 12 ) where U = e x e x + e y e y is the unit dyadic tensor 


in Eo. 


The contribution \d 1° Eq. (IBlbH can be expressed in terms of the pair correlation 
function p^ 2) ( 1,2) = p( 2 )(iq, an; r 2 , a 2 ) as 


-2tt 


c-2-k 


(B3) 


Xz)(ri 2 ) = 2tt/3p 2 / da L / da 2 p {2) ( l,2)sis 2 . 

Jo Jo 

For a homogeneous and isotropic fluid p( 2 )(l,2) = (p/27r) 2 p(l, 2). The normalized pair 
correlation function p(l,2) can then be expanded on the complete set of 2D orthogonal 
rotational invariants as in Refs 0.Q 


ft(l, 2) = 9 (1,2) - 1 = Y/ 2) , 


m,n 


where 


<h mn (l, 2) = cos(m/3i + n(3 2 ) . 


(B4) 


(B5) 


Pi = oi — 9 12 and /3 2 = a 2 — 0\ 2 are the angles of the vectors Si and Si and the direction 
r 12 4, 5] and the coefficients h m n (r 12 ) depends only on the sole modulus r 12 =|| r 12 ||. The 
two invariants D( 1, 2) = $ 11 ( 1 , 2) (minus the angular part of the dipole-dipole energy in the 
plane E 2 ) and A(l,2) = ^> 1 _ 1 (1,2) (scalar product Si • s 2 ) play an important role since we 
have, as a short calculation will show 

2 ,.2 

2 


Xd 


( r i 2 ) = ^ [h A ( i 2 )U + h D (r i 2 ) (2r'i 2 r'i 2 - U)] 
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where r ]2 = ri 2 /ri 2 . Gathering the intermediate results we obtain 


XO12) = 2/I( r i 2 ) + yh A (ri 2 )U + y/i D (r i2 ) (2ri 2 ri 2 - U) . (B6) 


We turn now our attention to the dipolar Green’s functions in E 2 . The bare Green’s 
function is given by 


r* , \ 1 d & 

G 0 (ri,r 2 ) = -^-^-\ogr l2 , 

27t or i ot 2 

— 7TZ~ (2 r 12 r 12 — U) . 


27r r 


(B7) 


12 


It is important to remark that Eq. (IB7D implies that 


Tr G 0 (ri, r 2 ) = ^-A(-logri 2 ) = -5 (2) ( r i2) • (B8) 

Z7T 

We note that in Fourier space Go = —kk (with k = k /k) identifies with minus the 
projector P/ = kk on longitudinal modes. Denoting by Q/, = U — P/, the projector on 
tranverse modes one has the relations P^ ■ P^ = P*,, Qfc • Q& = Q k and P k ■ Qfc = 0 from 
which the identity (aP^ + bQ k ) L = a _1 Pfc + b~ l Q k can be deduced. Assuming the locality 
of dielectric tensor is local, i.e. e — eU in Fourier space, permits to compute the dressed 
Green’s function G. One gets, as expected, 


G = -P fc [(1 + a) P k + cxQ,]- 1 = ^ (B9) 

In order to obtain an expression for the dielectric constant we take the trace of both 
sides of Fulton’s relation (174]) and integrate over the whole plane. Making use of the expres¬ 
sions (1B6D for x and the formula (1B8D of the trace of Go one finds 

(6 ~^ e + 1) = y(l + ^ f™ 2nrh A (r)dr ) (BIO) 

Fulton’s relation (174]) gives asymptotically (i.e. for r 12 >> £, range of the dielectric 
tensor) 

X( r i 2 ) ^ ^ g 1 ' ) 2 ^2 (2ri 2 ri 2 - U) . (Bll) 

A comparison of the above expression with Eq. (168]) yields readily the asymptotic behaviors 
of h D (r) and h^(r : 


h^-(r) 

h7 m (r) 


0, 

(e-1) 2 1 


1 

2 ‘ 
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(B12a) 

(B12b) 








from which we conclude that, while /?^ ym '(r) is a short range function, the projection 
/i)y m '(r) exhibits a long range algebraic asymptotic decay. The results obtained in this 
appendix were obtained by other methods in Refs. 0,1 


Appendix C: Dipole displacements on the sphere S 2 


The algorithm devised below adapts to the 2 D case the algorithm used in Ref. 32] for 
the 3 D DHS fluid on the hypersphere £3. The initial configuration of the Markov chain 
is obtained by sampling N vectors z* uniformly on the sphere 1S2. The initial positions of 
the point dipoles are thus OM ! = Rz l . One uses the spherical coordinates in the base 
(e x , e^, e 2 ), i.e. z* = (sin 9 l cos </?*, sin 9 l sin 93 *, cos9 t ) T with cost/?* = 2£j — 1 and <p l = 
27t^ 2 where £j and ££ are 21V random numbers G (0,1). It is convenient to complete the 
orthogonal local basis of spherical coordinates at point M l by defining also u* = dz 1 /d6 l = 
(cos#* cos</3*, cos#* sin— sin6 I *) T and v* = dz l /dip 1 / sin #* = (— sin cos ip 1 , 0) T . The 
initial dipole [i l = /is 1 is randomly sampled in the plane 7~(M*) tangent to the sphere at 
point M l according to s* = cos 0*u* + sin 0*v* with 0* = 2712(3 where £3 is a random number 

e(0 ,l). 

Due to the curvature of the space the trial move of dipole /T is made in three steps 
in which a displacement and a rotation are involved. First, the new position z„ ew of the 
(randomly chosen) dipole “i” is chosen uniformly on a small cap of angle <5# about point M* 
(i.e. the intersection of a 3D cone of axis z* and angle 56 with the sphere). It is convenient 
to use spherical coordinates in the local basis (z*, u*, v*) so as to define 

z’new = sin 9 l new cos <p l new u l + sin 6' new sin < ew v* + cos 9' new z* (Cl) 

The choice cos #( iew = (1 — cos 56)^1 + cos 59 and (p l new = 27r^ where £4 and £5 are random 
numbers G (0,1) ensures a uniform sampling of the cap. The new local basis at point z l new 
is then given by 


<ew = cos 9 l new cos </4 w u* + cos #^ ew sin </4w v i - sin 9* new z‘ 


(C2) 


V new = - sin (fi l new U* + COS <^ ew V* . (C3) 

The second step consits in rotating the vector s* in the tangent plane T(M») by an 
incremental angle <50* so that the new vector Sg tepl reads as 

(C4) 


s ste P i = cos <50*s* + sin <50jZ* x s* 
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where the choice 50j = (£g — 0.5)50 with £g some random number G (0,1) ensures a uniform 
sampling in the interval (—50/2,50/2). 

However vector Sg tepl does not belong to the plane T(M* ew ) and the third and last step 
of the trial MC move consists in a parallel transport of Sg tepl from the point M ! to the new 
a priori position M/ ew according to Eq. (EDD s/ ew = Tww r ; ew s;0 epl , i.e. 


^new ^Stepl 


^Stepl ^new 

1 + COS0* 


z 2 + z 


new) 


(C5) 


The trial displacement of dipole s* is then accepted or rejected according to the standard 
Metropolis criterion 33] after testing the overlaps and the energies. If the trial move is 
accepted we make the substitutions z] iew — > z\ uj, ew —* iT, v* ew —* v*, and s] iew — > s*. In the 
numerical experiments reported in Sec. |V] the parameters 50 and 50 were chosen in such a 
way that the acceptance ratio was ~ 0.5. 
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N 

/3 U mono 

^mono 

^00 mono/^j 

h A mono (cj) 

hP mono (cj) 

^■mono 

250 

-1.78533(8) 

18.14(3) 

4.335(1) 

2.821(2) 

4.580(2) 

3.954(2) 

500 

-1.78765(7) 

18.09(3) 

4.340(1) 

2.819(2) 

4.582(2) 

3.970(2) 

1000 

-1.78876(8) 

17.98(5) 

4.341(1) 

2.823(2) 

4.577(2) 

3.977(2) 

2000 

-1.78946(8) 

17.93(7) 

4.342(1) 

2.817(2) 

4.575(2) 

3.981(2) 

4000 

-1.78969(8) 

18.22(10) 

4.339(1) 

2.816(2) 

4.571(2) 

3.980(2) 

oo 

-1.78999(5) 

- 

- 

- 

- 

- 

N 

Pu hi 

e bi 

g00 U {(J ) 

h A bi (a) 

h D bi (cr) 

Z bi 

250 

-1.79010(7) 

18.06(2) 

4.334(1) 

2.800(2) 

4.583(2) 

3.962(2) 

500 

-1.78998(8) 

18.01(3) 

4.339(1) 

2.811(2) 

4.579(2) 

3.974(2) 

1000 

-1.79003(8) 

17.83(5) 

4.339(1) 

2.814(2) 

4.575(2) 

3.974(2) 

2000 

-1.79005(8) 

17.75(7) 

4.338(1) 

2.814(2) 

4.569(2) 

3.978(2) 

4000 

-1.79000(8) 

17.96(10) 

4.334(1) 

2.809(2) 

4.562(2) 

3.975(2) 

oo 

-1.79000(6) 

- 

- 

- 

- 

- 


TABLE I: Number of particles N, reduced internal energy per particle /3u, dielectric constant e, 
contact values g 00 {a ), h A (a), and h D (a) of some projections of the pair correlation function, and 
compressibility factor Z = (3P/p of the DHS fluid in the state (p* = 0.7, p* 2 = 2) for mono-dipoles 
(top) and bi-dipoles (bottom). The number in bracket corresponds to two standard deviations and 
gives the accuracy of the last digits. The thermodynamic limit of the internal energies were obtained 
from a linear regression in the variable N _1 . For each state lVc on f. = 16.0 x 10 9 configurations 
were generated. 
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N 

/3u mono 

^mono 

^00 mono( (J ) 

h A mono (a) 

h D mono( a ) 

^mono 

250 

-4.1625(3) 

43.94(9) 

5.143(2) 

5.531(4) 

7.146(3) 

1.660(2) 

500 

-4.1670(3) 

44.36(13) 

5.142(2) 

5.523(3) 

7.135(3) 

1.667(2) 

1000 

-4.1693(3) 

42.99(18) 

5.136(2) 

5.513(4) 

7.123(3) 

1.665(2) 

2000 

-4.1703(3) 

42.57(26) 

5.133(2) 

5.507(4) 

7.117(3) 

1.664(2) 

4000 

-4.1710(3) 

41.06(36) 

5.123(2) 

5.495(4) 

7.099(3) 

1.656(2) 

oo 

-4.17148(16) 

- 

- 

- 

- 

- 

N 

/3u bi 

e bi 

g°° bi (cr) 

h A bi (cr) 

h D hi (a) 

Z bi 

250 

-4.1722(3) 

44.55(9) 

5.136(2) 

5.501(3) 

7.141(3) 

1.657(2) 

500 

4.1718(3) 

43.96(11) 

5.136(2) 

5.506(4) 

7.131(3) 

1.657(2) 

1000 

-4.1717(3) 

43.24(18) 

5.132(2) 

5.502(4) 

7.117(3) 

1.662(2) 

2000 

-4.1714(3) 

42.89(27) 

5.122(2) 

5.487(4) 

7.103(3) 

1.655(2) 

4000 

-4.1714(3) 

45.14() 

5.137(2) 

5.510(4) 

7.123(3) 

1.670(2) 

oo 

-4.17138(15) 

- 

- 

- 

- 

- 


TABLE II: Same as Tabled but for the state (p* = 0.6, p* 2 = 4). 
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* *2 

(p = 0.7, (i = 2) 



* *2 

(p = 0.6, ji = 4) 



FIG. 1: Dimensionless MC internal energy /3u versus the inverse number of particles 1/N for the 
states ( p * = 0.7, ^ i* 2 = 2) (top) and ( p* = 0.6, /r* 2 = 4) (bottom). Black circles : mono-dipoles, 
red circles : bi-dipoles. The error bars correspond to two standard deviations. Blue dashed lines : 
linear regressions. 
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MONO-DIPOLES : (p*= 0.7, li“ 2 = 2) MONO-DIPOLES : (p = 0.6, Li 2 = 4) 




BI-DIPOLES : (p*= 0.7, p* 2 = 2) 


BI-DIPOLES : (p*= 0.6, p* 2 = 4) 



N 


conf 



FIG. 2: Cumulated dielectric constant e for the states I = (p* = 0.7, ^ i* 2 = 2) (left) and II 
= (p* = 0.6, /r* 2 = 4) (right) as a function of the number of configurations iV cori f. Top : mono¬ 


dipoles, bottom : bi-dipoles. 
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MONO-DIPOLES : (p*= 0.6, p* 2 = 4) 


MONO-DIPOLES : (p*= 0.7, a 2 = 2) 




BI-DIPOLES : (p*= 0.7, p* 2 = 2) BI-DIPOLES : (p*= 0.6, p >2 = 4) 




FIG. 3: Kirkwood’s factor G K for the states I=(p* = 0.7, //* 2 = 2) (left) and II=(/ 0 * = 0.6, /r* 2 = 4) 
(right) and various numbers N of particles as a function of the reduced distance r* = r/a. Top 
: mono-dipoles, bottom : bi-dipoles. Solid lines : MC data, dashed lines : asymptotic behaviors 
given by Eq. d78|) for nrono-dipoles and Eq. (IHHl) for bi-dipoles. 
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10 20 30 40 50 

j-* 


* *2 

(p = 0.6, 1 1 '= 4) 



FIG. 4: Projection h A (r*) for the states I=(/9* = 0.7, /i* 2 = 2) (top) and II=(/ 0 * = 0.6, fi* 2 = 4) 
(bottom) and various numbers N of particles (bi-dipoles) and N = 1000 (mono-dipoles). Solid 
lines : MC data, dashed lines : asymptotic behavior. 
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MONO-DIPOLES : (p = 0.7, p = 2) 


*2 


MONO-DIPOLES : (p = 0.6, p = 4) 




*2 


BI-DIPOLES : (p = 0.7, p = 2) 


*2 


BI-DIPOLES : (p = 0.6, p = 4) 




FIG. 5: Ratio h v (r*)/h^ symp (r*) for the state I=(p* = 0.7, p* 2 = 2) (left) and II=(/9* = 0.6, 
p* 2 = 4) (right) and various numbers N of particles. Top : mono-dipoles, bottom : bi-dipoles. 
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